library(dplyr)
library(ggplot2)
library(plotly)
library(glmnet)
library(gbm)
load("../1_data_management_dplyr/fichiers_prepared.RData")
nb_ET_CODGEO=finess_et%>%
group_by(CODGEO)%>%
summarise(nb_ET=n())
sum(!nb_ET_CODGEO$CODGEO%in%insee$CODGEO)
## [1] 3191
nb_ET_CODGEO=merge(nb_ET_CODGEO,insee,by="CODGEO",all.y=T)
nb_ET_CODGEO <- nb_ET_CODGEO%>%mutate(nb_ET=ifelse(is.na(nb_ET),0,nb_ET))
table(nb_ET_CODGEO$nb_ET)%>%head
##
## 0 1 2 3 4 5
## 28153 2888 1400 947 583 425
data_model <- nb_ET_CODGEO%>%
select(-LIBGEO,-CODGEO)%>%
mutate_if(is.character,as.numeric)
train=sample(x = 1:nrow(data_model),size = round(.65*nrow(data_model)))
Les valeurs manquantes sont très mal gérées par les GLM, on va donc imputer par la moyenne.
Il existe de nombreuses stratégies d’imputation donc certaines s’appuient sur des modèles de machine learning avancé : - Moyenne - Médiane - HotDeck - HotDeck stratifié - k plus proches voisins (kNN) - Régression/classification GLM - Régression/classification CART - Régression/classification GBM/RandomForest
Il existe plusieurs packages R qui implémente des stratégies d’imputation comme MICE ou simputation.
Lorsque les contraintes et le contexte sont très spécifiques et demande une maîtrise importante, il ne faut pas hésiter à implémenter son modèle d’imputation.
Vous pouvez vous référer à la présentation (Séminaire EIG #2) sur l’imputation des non-réponses partielles pour l’enquête OC : https://gitlab.com/DREES_code/OSAM/enquete_oc
data_model_imputed_for_glm=data_model%>%
mutate_all(function(x)ifelse(is.na(x),mean(x,na.rm=T),x))
Un GLM dans R ce décrit comme une formule Y~X1+X2+X3 où Y est la variable à expliquer/prédire et X1, X2, X4 les variables explicatives.
Lorsqu’on a déjà sélectionné les variables pertinentes pour le modèle, il suffit d’écrit . pour utiliser toutes les variables disponibles sauf Y.
Si on a souhaité conserver la variable ID et qu’on souhaite garder toutes les variables sauf ID il suffit d’écrire Y~.-ID.
On peut ajuster le modèle avec un certain nombre de parmaètres facultatifs. - Choisir la fonction de lien (log, inverse, logit, identité…) et la loi conditionnelle. Si on a oublié quelle est la fonction de lien canonique associée à une loi de la famille exponentielle, pas de problème, elle est assignée par défaut. - On peut préciser un ou plusieurs offset (coeff fixé, pas forcément à 1). - Pondérer les observations avec weights (par exemple si on veut donner moins d’importance aux données plus anciennes). Des idées générales sur cette approche, valable pour la plupart des modèles prédictifs. ici. - On peut définir les contrasts pour préciser quelle catégorie doit être utilisée comme référence pour une variable explicative catégorielle.
model <- glm(nb_ET~.,
data=data_model_imputed_for_glm[train,],
family=poisson(link = "log"))
coeff_glm=summary(model)$coefficients
coeff_glm
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.126276e+01 1.052242e+00 -10.703581 9.791915e-27
## NBMENFISC14 -5.527830e-05 1.609641e-06 -34.342000 1.854201e-258
## NBPERSMENFISC14 2.694129e-05 7.808639e-07 34.501902 7.511508e-261
## MED14 3.515499e-05 2.674837e-06 13.142852 1.870430e-39
## PIMP14 1.081355e-01 3.319893e-03 32.571972 1.023125e-232
## TP6014 -3.693168e-01 6.034702e-03 -61.198850 0.000000e+00
## TP60AGE114 -5.709656e-03 1.773697e-03 -3.219070 1.286070e-03
## TP60AGE214 -5.125060e-02 2.493470e-03 -20.553926 7.098366e-94
## TP60AGE314 -6.478064e-03 2.846824e-03 -2.275541 2.287349e-02
## TP60AGE414 1.785775e-02 3.142301e-03 5.683017 1.323390e-08
## TP60AGE514 2.139041e-02 3.350546e-03 6.384157 1.723447e-10
## TP60AGE614 -1.435137e-01 2.827824e-03 -50.750585 0.000000e+00
## TP60TOL114 -1.857351e-02 4.571160e-03 -4.063195 4.840554e-05
## TP60TOL214 -5.768398e-02 2.596958e-03 -22.212134 2.621841e-109
## PRA14 1.862361e-01 7.775149e-03 23.952742 8.651937e-127
## PTSAC14 -1.303710e-02 5.698780e-03 -2.287701 2.215494e-02
## PPEN14 2.039243e-01 9.969454e-03 20.454914 5.431799e-93
## PPAT14 2.629440e-01 9.632152e-03 27.298574 4.410173e-164
## PPSOC14 -2.113385e-01 1.078465e-01 -1.959624 5.003979e-02
## PPFAM14 7.672117e-01 1.076224e-01 7.128736 1.012950e-12
## PPMINI14 5.105354e-01 1.095562e-01 4.660033 3.161587e-06
## PPLOGT14 3.224468e+00 1.121298e-01 28.756559 7.499055e-182
## D114 -1.010405e-03 2.484651e-05 -40.665886 0.000000e+00
## D914 -1.485463e-04 5.893037e-06 -25.207092 3.349085e-140
## RD14 1.362044e+00 6.057189e-02 22.486407 5.638491e-112
# install.packages("stargazer")
stargazer::stargazer(model, type = "html")
##
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>nb_ET</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NBMENFISC14</td><td>-0.0001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">NBPERSMENFISC14</td><td>0.00003<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">MED14</td><td>0.00004<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PIMP14</td><td>0.108<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP6014</td><td>-0.369<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.006)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE114</td><td>-0.006<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE214</td><td>-0.051<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE314</td><td>-0.006<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE414</td><td>0.018<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE514</td><td>0.021<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60AGE614</td><td>-0.144<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60TOL114</td><td>-0.019<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.005)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">TP60TOL214</td><td>-0.058<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.003)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PRA14</td><td>0.186<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.008)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PTSAC14</td><td>-0.013<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.006)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PBEN14</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPEN14</td><td>0.204<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.010)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPAT14</td><td>0.263<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.010)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPSOC14</td><td>-0.211<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.108)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPFAM14</td><td>0.767<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.108)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPMINI14</td><td>0.511<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.110)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PPLOGT14</td><td>3.224<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.112)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">PIMPOT14</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">D114</td><td>-0.001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00002)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">D914</td><td>-0.0001<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">RD14</td><td>1.362<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.061)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-11.263<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.052)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>23,804</td></tr>
## <tr><td style="text-align:left">Log Likelihood</td><td>-60,696.850</td></tr>
## <tr><td style="text-align:left">Akaike Inf. Crit.</td><td>121,443.700</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
pred_glmtrain=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[train,],type="response"), obs=data_model[train,]$nb_ET)
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading
perf_train=data.frame(pred=data_model[train,]$nb_ET, obs=data_model[train,]$nb_ET)
pred_glmtest=data.frame(pred=predict(model,newdata = data_model_imputed_for_glm[-train,],type="response"),obs=data_model[-train,]$nb_ET)
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading
perf_test=data.frame(pred=data_model[-train,]$nb_ET, obs=data_model[-train,]$nb_ET)
Il s’agit des courbes de mesure des inégalités de richesse.
Ici on les applique à la mesure des inégalités de Y (dans notre cas nb_ET) dans la population : - La courbe “perfection” décrit les vraies d’inégalités - La courbe “random” décrit un modèle incapable de discerner les inégalités - La courbe “glm” décrit la capacité de notre modèle à appréhender (car échantillon de test) les inégalités
L’intérêt d’une telle mesure est qu’on peut définir - ce qu’est un “mauvais” modèle (random) - ce qu’est un “bon” modèle (perfection) et comparer notre modèle à ces deux extrêmes.
Cette approche est inspirée de l’étude de la courbe de ROC pour les modèles binaires (Bernoulli).
L’inconvénient est que cette métrique mesure les inégalités en rang et pas en taille.
pred_glmtest%>%
na.omit%>%
arrange(-pred)%>%
mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
{
data.frame(y100=quantile(.$y,0:100/100),
x100=quantile(.$x,0:100/100))
}->Lorenz_glm_points
pred_glmtest%>%
na.omit%>%
arrange(-obs)%>%
mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
{
data.frame(y100=quantile(.$y,0:100/100),
x100=quantile(.$x,0:100/100))
}->Perfection
g <- ggplot()+
geom_line(data = Lorenz_glm_points,
aes(x=x100,y=y100,color="glm"))+
geom_line(data = Perfection,
aes(x=x100,y=y100,color="perfection")) +
geom_line(data=data.frame(x100=c(0,1), yrandom=c(0,1)),
aes(x=x100,y=yrandom,color="random"))
g
Par extension de l’AUC sous la ROC, on calcule l’AUC sous la courbe de Lorenz.
Cette métrique, translatée (X->2*X-1) pour se comparer à “random” est appelée indice de Gini.
Pour bien faire, on peut rapporter le Gini du modèle au Gini de la “perfection” pour que notre indice ait des valeurs POSSIBLES entre 0 et 1.
Gini=function(pred){
pred%>%
na.omit%>%
arrange(-pred)%>%
mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
{
# print(paste0("nombre_de_lignes:",nrow(.)," indice de Gini:"))
2*mean(.$y)-1
}
}
Gini(pred_glmtest)
## [1] 0.5972249
Gini(perf_test)
## [1] 0.9289132
Gini(pred_glmtrain)
## [1] 0.6228157
Gini(perf_train)
## [1] 0.9325213
Performance faible mais pas d’overfitting ie peu d’écart entre apprentissage et test.
Très bonne explication ici. 2 paramètres à choisir : - alpha qui équilibre entre les pénalisations L1 et L2. - lambda qui équilibre la fonction de perte entre vraisemblance et pénalisation sur la taille des coefficients Lambda élevé signifie une pénalisation importante sur la taille des coefficients.
Lambda faible signifie qu’on est proche d’un GLM classique non pénalisé.
Le package R glmnet test automatiquement et de manière optimisé plusieurs lambda.
En revanche c’est à nous de tester plusieurs paramètres alpha. C’est ce qu’on appelle l’optimisation des hyperparamètres (hyperparameters tuning). La méthode qui consiste à tester les combinaisons des différents hyperparamètres s’appelle gridsearch.
La méthode gourmande/gloutonne (bruteforce) consiste à tester toutes les combinaisons.
Il existe des méthodes bayesiennes pour explorer intelligemment la grille de valeurs.
De nombreux packages R sont disponibles, mlrMBO est plutôt bien maintenu.
Mais ces considérations dépassent l’ambition de cette formation.
Dans R il y a beaucoup de liberté quant au traitement des données avec les formats matrix, data.frame, data.table…
Pour glm, on avait un modèle de la forme cible ~ variables explicatives ou Y~X avec Y et X dans un même data.frame.
Pour glmnet, on doit fournir Y comme un vecteur et X comme une matrice de valeurs numériques.
Par conséquent les variables catégorielles devront être binarisées (dummy variables), ce procédé s’appelle le one-hot encoding. C’est très simple en R, il suffit d’utiliser la fonction model.matrix. Un exemple ici
Pas besoin d’y recourir dans notre cas. Si
target_train = data_model_imputed_for_glm[train,]$nb_ET
explanatory_train = data_model_imputed_for_glm[train,]%>%
select(-nb_ET)%>%as.matrix
target_test = data_model_imputed_for_glm[-train,]$nb_ET
explanatory_test = data_model_imputed_for_glm[-train,]%>%
select(-nb_ET)%>%as.matrix
thresh est le coefficient de convergence. On peut lire dans la documentation :
thresh Convergence threshold for coordinate descent. Each inner coordinate-descent
loop continues until the maximum change in the objective after any coefficient
update is less than thresh times the null deviance. Defaults value is 1E-7.
glmnet_model <- glmnet(x=explanatory_train,
y=target_train,
family="poisson",
alpha=1,
nlambda=100,#par défaut
thresh = 1e-06,#par défaut
maxit=1e7)
plot(glmnet_model$lambda,glmnet_model$dev.ratio)
# Distribution des coefficients en fonction du lambda.
plot(glmnet_model, xvar='lambda',label=T)
s=sample(glmnet_model$lambda,1)
smin=min(glmnet_model$lambda)
pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
print(paste("s=0",Gini(pred_glmnettest)))
## [1] "s=0 0.591052172672043"
pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
print(paste("s=smin=",smin,"Gini=",Gini(pred_glmnettest)))
## [1] "s=smin= 0.00156900416023852 Gini= 0.591052172672043"
pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=s)[,1],obs=target_test)
print(paste("s=random",s,"Gini=",Gini(pred_glmnettest)))
## [1] "s=random 0.00249829832594608 Gini= 0.590321105877747"
Pour automatiser et rationnaliser le choix du lambda on fait de la validation croisée (cross-validation) ie on coupe le dataset en NFOLDS, disons en 10, et entraîne sur 90% vs validation sur 10% 10 fois.
NFOLDS=10
for (alpha in 0:10/10){
print(alpha)
tryCatch(rm(glmnet_cv,pred_glmnettest),warning = function(w) {
print(paste("warning",w))
}, error = function(e) {
print(paste("error",e))
})
tryCatch({
glmnet_cv = cv.glmnet(x = explanatory_train, y = target_train,
family = "poisson",#'gaussian',
# Pénalisation L1 100%
alpha = alpha,
# On s'intéresse à la deviance - on suppose que la distribution conditionnelle suit une loi de Poisson
type.measure = "deviance",#'rmse'
# 5-fold cross-validation
nfolds = NFOLDS,
# valeurs élevée pour un entraînement plus rapide mais moins performant
thresh = 1e-3,
# On limite le nombre d'itérations
maxit = 1e5)
s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=0)[,1],obs=target_test)
print(paste("s=0",Gini(pred_glmnettest)))
pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
print(paste("s=smin",Gini(pred_glmnettest)))
pred_glmnettest=data.frame(pred=predict(glmnet_cv,newx = explanatory_test,type="response",s=s1se)[,1],obs=target_test)
print(paste("s=s1se",Gini(pred_glmnettest)))
},warning = function(w) {
print(paste("warning",w))
}, error = function(e) {
print(paste("error",e))
})
}
## [1] 0
## [1] "warning simpleWarning in rm(glmnet_cv, pred_glmnettest): objet 'glmnet_cv' introuvable\n"
## [1] "s=0 0.530549879254173"
## [1] "s=smin 0.491033845212646"
## [1] "s=s1se -0.125987883560982"
## [1] 0.1
## [1] "s=0 0.569958405168113"
## [1] "s=smin 0.567633591437849"
## [1] "s=s1se 0.543620369155664"
## [1] 0.2
## [1] "s=0 0.570660871328051"
## [1] "s=smin 0.56847105691447"
## [1] "s=s1se 0.556439987548278"
## [1] 0.3
## [1] "s=0 0.571030222945907"
## [1] "s=smin 0.569235536742053"
## [1] "s=s1se 0.547383521034902"
## [1] 0.4
## [1] "s=0 0.571227063029675"
## [1] "s=smin 0.569423306750592"
## [1] "s=s1se 0.546153643328189"
## [1] 0.5
## [1] "s=0 0.57165069821929"
## [1] "s=smin 0.569992293677941"
## [1] "s=s1se 0.543537715483627"
## [1] 0.6
## [1] "s=0 0.571429104538633"
## [1] "s=smin 0.569558630255826"
## [1] "s=s1se 0.594412732800699"
## [1] 0.7
## [1] "s=0 0.571385886602735"
## [1] "s=smin 0.569975162114087"
## [1] "s=s1se 0.546906481184709"
## [1] 0.8
## [1] "s=0 0.571488474268542"
## [1] "s=smin 0.569781563915826"
## [1] "s=s1se 0.551642407637421"
## [1] 0.9
## [1] "s=0 0.571581127356319"
## [1] "s=smin 0.56966485603759"
## [1] "s=s1se 0.550196447455412"
## [1] 1
## [1] "s=0 0.572030488708477"
## [1] "s=smin 0.570142709961255"
## [1] "s=s1se 0.548101173627391"
Gini(pred_glmtest)
## [1] 0.5972249
plot(glmnet_cv)
Des idées pratiques en anglais ici
et théoriques en français ici
s1se=glmnet_cv$lambda.1se
smin=glmnet_cv$lambda.min
coeff_glmnet=coef(glmnet_model,s=s1se)
coeff_glmnet
## 27 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 1.353915e+01
## NBMENFISC14 -4.239680e-08
## NBPERSMENFISC14 2.691878e-07
## MED14 1.551383e-05
## PIMP14 7.410487e-02
## TP6014 -2.098730e-01
## TP60AGE114 .
## TP60AGE214 -4.625761e-02
## TP60AGE314 -4.476880e-04
## TP60AGE414 4.092152e-03
## TP60AGE514 8.465594e-03
## TP60AGE614 -1.363710e-01
## TP60TOL114 -4.754384e-02
## TP60TOL214 -7.678858e-02
## PRA14 -1.093402e-02
## PTSAC14 .
## PBEN14 .
## PPEN14 .
## PPAT14 2.277342e-02
## PPSOC14 .
## PPFAM14 7.653354e-02
## PPMINI14 .
## PPLOGT14 2.198343e+00
## PIMPOT14 -1.571213e-01
## D114 -1.243651e-03
## D914 -1.234753e-05
## RD14 .
# On prépare la table des coeffs de glmnet
coeff_glmnet=as.matrix(coeff_glmnet)
coeff_glmnet=data.frame(var=rownames(coeff_glmnet),coeff_glmnet=coeff_glmnet[,1])
coeff_glmnet$var=as.character(coeff_glmnet$var)
# On prépare la table des coeffs de glm
coeff_glm=data.frame(var=rownames(coeff_glm),coeff_glm=coeff_glm[,1])
coeff_glm$var=as.character(coeff_glm$var)
# On joint les deux tables et on les compare
comparaison_coeff=merge(coeff_glm,coeff_glmnet,by="var")
comparaison_coeff%>%
mutate(ratio_glm_sur_glmnet=coeff_glm/coeff_glmnet)%>%
as.tbl
## # A tibble: 25 x 4
## var coeff_glm coeff_glmnet ratio_glm_sur_glmnet
## <chr> <dbl> <dbl> <dbl>
## 1 (Intercept) -11.3 1.35e+1 -0.832
## 2 D114 -0.00101 -1.24e-3 0.812
## 3 D914 -0.000149 -1.23e-5 12.0
## 4 MED14 0.0000352 1.55e-5 2.27
## 5 NBMENFISC14 -0.0000553 -4.24e-8 1304.
## 6 NBPERSMENFISC14 0.0000269 2.69e-7 100.
## 7 PIMP14 0.108 7.41e-2 1.46
## 8 PPAT14 0.263 2.28e-2 11.5
## 9 PPEN14 0.204 0. Inf
## 10 PPFAM14 0.767 7.65e-2 10.0
## # ... with 15 more rows
Pour bien comprendre ce qui se passe, on jette un oeil à la matrice des corrélations. On se rend compte que lorsque deux variables sont très corrélées, en général le glmnet n’apprend que sur l’une des deux ie le coeff de l’autre passe à 0 ! C’est de la sélection de variables
cor(explanatory_train)%>%as.data.frame->cormat
rownames(cormat) <- colnames(cormat)
knitr::kable(cormat)
| NBMENFISC14 | NBPERSMENFISC14 | MED14 | PIMP14 | TP6014 | TP60AGE114 | TP60AGE214 | TP60AGE314 | TP60AGE414 | TP60AGE514 | TP60AGE614 | TP60TOL114 | TP60TOL214 | PRA14 | PTSAC14 | PBEN14 | PPEN14 | PPAT14 | PPSOC14 | PPFAM14 | PPMINI14 | PPLOGT14 | PIMPOT14 | D114 | D914 | RD14 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NBMENFISC14 | 1.0000000 | 0.9980597 | 0.0215953 | 0.0028921 | 0.1052056 | -0.0860245 | -0.0521292 | 0.0070649 | 0.0158168 | -0.0142194 | -0.1186299 | -0.0522551 | 0.0150650 | 0.0525470 | 0.0455294 | 0.0217454 | -0.0628268 | 0.0333159 | 0.0691816 | -0.0199791 | 0.0866507 | 0.1057856 | -0.1308763 | -0.1385921 | 0.1068657 | 0.2617827 |
| NBPERSMENFISC14 | 0.9980597 | 1.0000000 | 0.0229011 | 0.0019599 | 0.1172289 | -0.0808335 | -0.0463743 | 0.0150583 | 0.0273077 | 0.0018635 | -0.1046563 | -0.0389131 | 0.0225112 | 0.0648001 | 0.0592228 | 0.0118646 | -0.0786364 | 0.0229949 | 0.0823872 | -0.0010762 | 0.0958551 | 0.1146606 | -0.1300055 | -0.1458092 | 0.1049247 | 0.2665857 |
| MED14 | 0.0215953 | 0.0229011 | 1.0000000 | 0.4077471 | -0.2910832 | -0.1500026 | -0.1992152 | -0.2260350 | -0.1888196 | -0.1312729 | -0.0769175 | -0.1583395 | -0.2742843 | 0.2530778 | 0.2292955 | 0.0560573 | -0.1551675 | 0.1621254 | -0.3739791 | -0.3134540 | -0.3341195 | -0.3606243 | -0.3806792 | 0.3803337 | 0.4262154 | 0.1411290 |
| PIMP14 | 0.0028921 | 0.0019599 | 0.4077471 | 1.0000000 | -0.6622976 | -0.3469204 | -0.4525352 | -0.5028118 | -0.4179267 | -0.2907931 | -0.1771124 | -0.3892491 | -0.6180156 | 0.5867046 | 0.5423128 | 0.0777533 | -0.2849247 | 0.1704712 | -0.8006538 | -0.5981922 | -0.7543827 | -0.7850175 | -0.8155333 | 0.8272099 | 0.6691432 | 0.0472803 |
| TP6014 | 0.1052056 | 0.1172289 | -0.2910832 | -0.6622976 | 1.0000000 | 0.4284438 | 0.5973560 | 0.7048388 | 0.5996171 | 0.4650306 | 0.2822975 | 0.5854329 | 0.8107210 | -0.3231671 | -0.2910016 | -0.0803133 | 0.0956991 | -0.2073212 | 0.8369567 | 0.5302835 | 0.8542956 | 0.8167357 | 0.4226966 | -0.7394627 | -0.3623792 | 0.2468963 |
| TP60AGE114 | -0.0860245 | -0.0808335 | -0.1500026 | -0.3469204 | 0.4284438 | 1.0000000 | 0.6996132 | 0.5802272 | 0.6210235 | 0.4895079 | 0.3840939 | 0.4727018 | 0.4755313 | -0.2137757 | -0.2080100 | 0.0222543 | 0.0734345 | -0.1167832 | 0.4565836 | 0.3471938 | 0.4583764 | 0.3999994 | 0.3083101 | -0.2192054 | -0.2596827 | -0.0823605 |
| TP60AGE214 | -0.0521292 | -0.0463743 | -0.1992152 | -0.4525352 | 0.5973560 | 0.6996132 | 1.0000000 | 0.8224020 | 0.7893669 | 0.5453349 | 0.3082232 | 0.5997015 | 0.6355497 | -0.2753687 | -0.2578959 | -0.0201535 | 0.1056694 | -0.1578076 | 0.5852137 | 0.4322669 | 0.5901522 | 0.5218786 | 0.3844230 | -0.3368280 | -0.3315089 | -0.0611591 |
| TP60AGE314 | 0.0070649 | 0.0150583 | -0.2260350 | -0.5028118 | 0.7048388 | 0.5802272 | 0.8224020 | 1.0000000 | 0.7907886 | 0.5511260 | 0.3137646 | 0.6514217 | 0.7346698 | -0.2727765 | -0.2521532 | -0.0360734 | 0.0860176 | -0.1825043 | 0.6461682 | 0.4656216 | 0.6521336 | 0.5870576 | 0.4062684 | -0.4339575 | -0.3475086 | 0.0125118 |
| TP60AGE414 | 0.0158168 | 0.0273077 | -0.1888196 | -0.4179267 | 0.5996171 | 0.6210235 | 0.7893669 | 0.7907886 | 1.0000000 | 0.7307591 | 0.4423994 | 0.7224931 | 0.5913820 | -0.1538919 | -0.1390782 | -0.0357988 | -0.0263840 | -0.1879764 | 0.5598322 | 0.4517782 | 0.5675372 | 0.4582539 | 0.3626659 | -0.3388408 | -0.2962708 | 0.0043506 |
| TP60AGE514 | -0.0142194 | 0.0018635 | -0.1312729 | -0.2907931 | 0.4650306 | 0.4895079 | 0.5453349 | 0.5511260 | 0.7307591 | 1.0000000 | 0.6496638 | 0.7462134 | 0.4171335 | -0.0159521 | -0.0118301 | -0.0162803 | -0.1432056 | -0.1691843 | 0.4257486 | 0.3814942 | 0.4536359 | 0.2806997 | 0.2724857 | -0.2270215 | -0.1962145 | 0.0273695 |
| TP60AGE614 | -0.1186299 | -0.1046563 | -0.0769175 | -0.1771124 | 0.2822975 | 0.3840939 | 0.3082232 | 0.3137646 | 0.4423994 | 0.6496638 | 1.0000000 | 0.6533659 | 0.2502127 | 0.0232300 | 0.0131133 | 0.0437007 | -0.1285273 | -0.1326114 | 0.2720856 | 0.2327039 | 0.3504872 | 0.1039938 | 0.1872943 | -0.1129029 | -0.1122044 | 0.0167602 |
| TP60TOL114 | -0.0522551 | -0.0389131 | -0.1583395 | -0.3892491 | 0.5854329 | 0.4727018 | 0.5997015 | 0.6514217 | 0.7224931 | 0.7462134 | 0.6533659 | 1.0000000 | 0.5554042 | -0.1217988 | -0.1263670 | 0.0508431 | -0.0506555 | -0.1472988 | 0.4968401 | 0.3935647 | 0.5602870 | 0.3325648 | 0.3090969 | -0.3250373 | -0.2132905 | 0.0812710 |
| TP60TOL214 | 0.0150650 | 0.0225112 | -0.2742843 | -0.6180156 | 0.8107210 | 0.4755313 | 0.6355497 | 0.7346698 | 0.5913820 | 0.4171335 | 0.2502127 | 0.5554042 | 1.0000000 | -0.3389176 | -0.3123168 | -0.0495672 | 0.1307905 | -0.2078403 | 0.7202636 | 0.5322553 | 0.7134434 | 0.6614165 | 0.4911544 | -0.5715571 | -0.4330499 | 0.0210353 |
| PRA14 | 0.0525470 | 0.0648001 | 0.2530778 | 0.5867046 | -0.3231671 | -0.2137757 | -0.2753687 | -0.2727765 | -0.1538919 | -0.0159521 | 0.0232300 | -0.1217988 | -0.3389176 | 1.0000000 | 0.9789818 | -0.1330281 | -0.8648094 | -0.3531594 | -0.3285947 | 0.0254452 | -0.4175110 | -0.4277412 | -0.4454008 | 0.4861982 | 0.3883575 | 0.0518605 |
| PTSAC14 | 0.0455294 | 0.0592228 | 0.2292955 | 0.5423128 | -0.2910016 | -0.2080100 | -0.2578959 | -0.2521532 | -0.1390782 | -0.0118301 | 0.0131133 | -0.1263670 | -0.3123168 | 0.9789818 | 1.0000000 | -0.3323673 | -0.8615756 | -0.4084203 | -0.2677540 | 0.0874546 | -0.3730935 | -0.3650827 | -0.3557900 | 0.4478249 | 0.3207720 | 0.0099472 |
| PBEN14 | 0.0217454 | 0.0118646 | 0.0560573 | 0.0777533 | -0.0803133 | 0.0222543 | -0.0201535 | -0.0360734 | -0.0357988 | -0.0162803 | 0.0437007 | 0.0508431 | -0.0495672 | -0.1330281 | -0.3323673 | 1.0000000 | 0.1876612 | 0.3516002 | -0.2183883 | -0.3073275 | -0.1176672 | -0.2039064 | -0.3307316 | 0.0721418 | 0.2371115 | 0.1914874 |
| PPEN14 | -0.0628268 | -0.0786364 | -0.1551675 | -0.2849247 | 0.0956991 | 0.0734345 | 0.1056694 | 0.0860176 | -0.0263840 | -0.1432056 | -0.1285273 | -0.0506555 | 0.1307905 | -0.8648094 | -0.8615756 | 0.1876612 | 1.0000000 | 0.1584739 | 0.0784269 | -0.2469206 | 0.1906209 | 0.2033623 | 0.1828383 | -0.2706725 | -0.2480069 | -0.0825585 |
| PPAT14 | 0.0333159 | 0.0229949 | 0.1621254 | 0.1704712 | -0.2073212 | -0.1167832 | -0.1578076 | -0.1825043 | -0.1879764 | -0.1691843 | -0.1326114 | -0.1472988 | -0.2078403 | -0.3531594 | -0.4084203 | 0.3516002 | 0.1584739 | 1.0000000 | -0.3916229 | -0.5438205 | -0.2652198 | -0.2918097 | -0.3397489 | 0.1757270 | 0.4534294 | 0.3265344 |
| PPSOC14 | 0.0691816 | 0.0823872 | -0.3739791 | -0.8006538 | 0.8369567 | 0.4565836 | 0.5852137 | 0.6461682 | 0.5598322 | 0.4257486 | 0.2720856 | 0.4968401 | 0.7202636 | -0.3285947 | -0.2677540 | -0.2183883 | 0.0784269 | -0.3916229 | 1.0000000 | 0.7775157 | 0.9418766 | 0.9502929 | 0.6102688 | -0.8001881 | -0.5919658 | 0.0197915 |
| PPFAM14 | -0.0199791 | -0.0010762 | -0.3134540 | -0.5981922 | 0.5302835 | 0.3471938 | 0.4322669 | 0.4656216 | 0.4517782 | 0.3814942 | 0.2327039 | 0.3935647 | 0.5322553 | 0.0254452 | 0.0874546 | -0.3073275 | -0.2469206 | -0.5438205 | 0.7775157 | 1.0000000 | 0.5513702 | 0.6205474 | 0.6262319 | -0.4804263 | -0.6514635 | -0.3087437 |
| PPMINI14 | 0.0866507 | 0.0958551 | -0.3341195 | -0.7543827 | 0.8542956 | 0.4583764 | 0.5901522 | 0.6521336 | 0.5675372 | 0.4536359 | 0.3504872 | 0.5602870 | 0.7134434 | -0.4175110 | -0.3730935 | -0.1176672 | 0.1906209 | -0.2652198 | 0.9418766 | 0.5513702 | 1.0000000 | 0.9076897 | 0.5034243 | -0.7954744 | -0.4576716 | 0.1678671 |
| PPLOGT14 | 0.1057856 | 0.1146606 | -0.3606243 | -0.7850175 | 0.8167357 | 0.3999994 | 0.5218786 | 0.5870576 | 0.4582539 | 0.2806997 | 0.1039938 | 0.3325648 | 0.6614165 | -0.4277412 | -0.3650827 | -0.2039064 | 0.2033623 | -0.2918097 | 0.9502929 | 0.6205474 | 0.9076897 | 1.0000000 | 0.5383898 | -0.8379558 | -0.5243251 | 0.1148539 |
| PIMPOT14 | -0.1308763 | -0.1300055 | -0.3806792 | -0.8155333 | 0.4226966 | 0.3083101 | 0.3844230 | 0.4062684 | 0.3626659 | 0.2724857 | 0.1872943 | 0.3090969 | 0.4911544 | -0.4454008 | -0.3557900 | -0.3307316 | 0.1828383 | -0.3397489 | 0.6102688 | 0.6262319 | 0.5034243 | 0.5383898 | 1.0000000 | -0.5297794 | -0.8423480 | -0.4658049 |
| D114 | -0.1385921 | -0.1458092 | 0.3803337 | 0.8272099 | -0.7394627 | -0.2192054 | -0.3368280 | -0.4339575 | -0.3388408 | -0.2270215 | -0.1129029 | -0.3250373 | -0.5715571 | 0.4861982 | 0.4478249 | 0.0721418 | -0.2706725 | 0.1757270 | -0.8001881 | -0.4804263 | -0.7954744 | -0.8379558 | -0.5297794 | 1.0000000 | 0.5020017 | -0.2934796 |
| D914 | 0.1068657 | 0.1049247 | 0.4262154 | 0.6691432 | -0.3623792 | -0.2596827 | -0.3315089 | -0.3475086 | -0.2962708 | -0.1962145 | -0.1122044 | -0.2132905 | -0.4330499 | 0.3883575 | 0.3207720 | 0.2371115 | -0.2480069 | 0.4534294 | -0.5919658 | -0.6514635 | -0.4576716 | -0.5243251 | -0.8423480 | 0.5020017 | 1.0000000 | 0.6621550 |
| RD14 | 0.2617827 | 0.2665857 | 0.1411290 | 0.0472803 | 0.2468963 | -0.0823605 | -0.0611591 | 0.0125118 | 0.0043506 | 0.0273695 | 0.0167602 | 0.0812710 | 0.0210353 | 0.0518605 | 0.0099472 | 0.1914874 | -0.0825585 | 0.3265344 | 0.0197915 | -0.3087437 | 0.1678671 | 0.1148539 | -0.4658049 | -0.2934796 | 0.6621550 | 1.0000000 |
pred_glmnettrain=data.frame(pred=predict(glmnet_model,newx = explanatory_train,type="response")[,1],obs=target_train)
pred_glmnettest=data.frame(pred=predict(glmnet_model,newx = explanatory_test,type="response",s=smin)[,1],obs=target_test)
pred_glmnettest%>%
na.omit%>%
arrange(-pred)%>%
mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
{
data.frame(y100=quantile(.$y,0:100/100),
x100=quantile(.$x,0:100/100))
}->Lorenz_glmnet_points
g <- g +
geom_line(data = Lorenz_glmnet_points,
aes(x=x100,y=y100,color="glmnet"))
g%>%ggplotly
Gini(pred_glmtest)
## [1] 0.5972249
Gini(pred_glmnettest)
## [1] 0.5910522
Gini(perf_test)
## [1] 0.9289132
Gini(pred_glmtrain)
## [1] 0.6228157
Gini(pred_glmnettrain)
## [1] 0.0150469
Gini(perf_train)
## [1] 0.9325213
Commençons par utiliser un arbre de regression pour expliquer simplement le nombre d’établissements dans la commune. Ceci nous permettra de détecter les variables les plus influentes.
library(rpart)
options_rpart=rpart.control(minbucket =100, cp = 1E-3, maxdepth = 5)
tree <- rpart(nb_ET~., data=data_model_imputed_for_glm[train,],control = options_rpart)
rpart.plot::rpart.plot(tree,roundint=FALSE)
On offset la variable NBMENFISC14 pour calculer un nombre d’établissement par milliers de ménages. C’est aussi l’occasion d’introduire une autre visualisation plus interactive.
tree <- rpart(1000*nb_ET/NBMENFISC14~.-NBMENFISC14-NBPERSMENFISC14, data=data_model_imputed_for_glm[train,],control = options_rpart)
visNetwork::visTree(tree, main = "Arbre de décision interactif",
width = "80%", height = "400px")
L’avantage des modèles basés sur des arbres de décision, c’est qu’un arbre de décision sait bien gérer les NA.
Par exemple pour un arbre binaire, on peut imaginer séparer les NA dans un 3ème split, ou bien séparer NA vs le reste. La contrainte imposée par les variables continues est bien moins forte que dans un GLM. Avec un arbre si on coupe X à un seuil th (choisi par le modèle), on peut déciser de regrouper les NA avec le groupe 1 X<th ou avec le groupe 2 X>th.
L’imputation n’est pas nécessaire ici.
data_model%>%head
## nb_ET NBMENFISC14 NBPERSMENFISC14 MED14 PIMP14 TP6014 TP60AGE114
## 1 0 304 799.5 21576.7 NA NA NA
## 2 0 100 235.5 21672.9 NA NA NA
## 3 0 6087 13660.5 19756.1 56.8215 15.7534 19.4181
## 4 0 603 1661.5 23204.8 NA NA NA
## 5 0 48 102.0 22157.5 NA NA NA
## 6 0 1030 2635.0 21679.3 61.7476 8.8425 NA
## TP60AGE214 TP60AGE314 TP60AGE414 TP60AGE514 TP60AGE614 TP60TOL114
## 1 NA NA NA NA NA NA
## 2 NA NA NA NA NA NA
## 3 19.5204 19.1982 14.7159 NA NA 5.40116
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 NA NA NA NA NA NA
## TP60TOL214 PRA14 PTSAC14 PBEN14 PPEN14 PPAT14 PPSOC14 PPFAM14 PPMINI14
## 1 NA NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA
## 3 24.796 71.8 67.9 3.9 27.3 10.1 6.5 2.8 1.8
## 4 NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA
## 6 NA 76.3 73.8 2.5 26.5 8.2 4.3 2.7 0.8
## PPLOGT14 PIMPOT14 D114 D914 RD14
## 1 NA NA NA NA NA
## 2 NA NA NA NA NA
## 3 1.9 -15.7 10545.7 33376.5 3.16495
## 4 NA NA NA NA NA
## 5 NA NA NA NA NA
## 6 0.8 -15.3 12797.2 34818.9 2.72082
params=c(shrinkage=.01,
nb_trees=1000,
depth=10,
nminobs=10,
bag.frac=.2)
gbm_model <- gbm(nb_ET~.
,data=data_model[train,]
,verbose=T
,train.fraction=.8
,shrinkage = params[1]
,n.trees = params[2]
,interaction.depth = params[3]
,n.minobsinnode = params[4]
,bag.fraction = params[5]
)
## Distribution not specified, assuming gaussian ...
## Iter TrainDeviance ValidDeviance StepSize Improve
## 1 161.6112 124.9445 0.0100 1.1945
## 2 159.5060 122.9226 0.0100 1.8978
## 3 158.3081 121.8718 0.0100 1.2360
## 4 156.5109 120.2592 0.0100 1.7279
## 5 154.6768 118.8279 0.0100 1.6574
## 6 152.7859 117.3525 0.0100 1.5123
## 7 151.5436 116.1283 0.0100 1.2773
## 8 150.0042 114.8612 0.0100 1.4503
## 9 147.9612 113.2020 0.0100 1.6624
## 10 146.6657 112.0619 0.0100 1.3079
## 20 135.3937 103.1946 0.0100 1.2281
## 40 115.0836 86.8909 0.0100 0.8974
## 60 101.1656 76.1005 0.0100 0.4093
## 80 91.7919 70.0738 0.0100 0.4968
## 100 84.9559 66.0688 0.0100 0.2363
## 120 80.7105 63.8572 0.0100 0.0554
## 140 76.8298 62.0454 0.0100 0.1331
## 160 72.9588 60.7204 0.0100 -0.1268
## 180 70.0501 60.2225 0.0100 0.0747
## 200 68.0310 60.3504 0.0100 0.1119
## 220 66.0278 60.5001 0.0100 0.0799
## 240 64.1627 60.8294 0.0100 -0.0563
## 260 63.4552 60.9274 0.0100 -0.0511
## 280 62.1450 61.4469 0.0100 -0.0137
## 300 61.1630 61.8902 0.0100 -0.0064
## 320 60.1808 62.1462 0.0100 -0.0569
## 340 59.3843 62.2452 0.0100 0.0246
## 360 58.6128 62.4991 0.0100 -0.0155
## 380 57.9869 62.9084 0.0100 -0.0024
## 400 57.0480 63.4925 0.0100 0.0210
## 420 56.5894 63.7226 0.0100 0.0317
## 440 55.9468 63.9451 0.0100 -0.0220
## 460 55.5446 63.9862 0.0100 -0.0205
## 480 54.9475 64.5412 0.0100 -0.0437
## 500 54.6453 64.6200 0.0100 -0.0058
## 520 54.2332 65.5485 0.0100 -0.0352
## 540 53.5622 66.1040 0.0100 -0.0117
## 560 52.8856 66.4188 0.0100 0.0406
## 580 52.2394 66.3664 0.0100 0.0133
## 600 51.7129 66.8059 0.0100 -0.0268
## 620 51.3457 66.8829 0.0100 -0.0309
## 640 50.8133 67.4678 0.0100 -0.0568
## 660 50.4482 67.2470 0.0100 -0.0526
## 680 49.9401 67.5536 0.0100 0.0125
## 700 49.4829 67.7143 0.0100 -0.0252
## 720 49.0172 68.7258 0.0100 -0.0031
## 740 48.5996 69.1680 0.0100 -0.0201
## 760 48.2110 69.2517 0.0100 -0.0419
## 780 47.8247 69.2563 0.0100 -0.0649
## 800 47.4762 69.2609 0.0100 -0.1056
## 820 47.0627 69.5084 0.0100 -0.0395
## 840 46.6427 69.5816 0.0100 -0.0627
## 860 46.4269 69.1736 0.0100 -0.0253
## 880 46.0515 69.3142 0.0100 -0.0526
## 900 45.8251 69.4874 0.0100 -0.0566
## 920 45.1928 69.9097 0.0100 -0.0063
## 940 44.9353 70.4944 0.0100 -0.0054
## 960 44.7112 70.7408 0.0100 -0.0065
## 980 44.4438 70.9160 0.0100 0.0046
## 1000 44.1010 71.2983 0.0100 -0.0346
pred_gbmtrain=data.frame(pred=predict(gbm_model,
newdata = data_model[train,]),
obs=target_train)
## Using 190 trees...
pred_gbmtest=data.frame(pred=predict(gbm_model,
newdata = data_model[-train,]),
obs=target_test)
## Using 190 trees...
print(paste("gbm test:",Gini(pred_gbmtest)))
## [1] "gbm test: 0.846249444647002"
print(paste("gbm train:",Gini(pred_gbmtrain)))
## [1] "gbm train: 0.862142145815112"
print(paste("glm:",Gini(pred_glmnettest)))
## [1] "glm: 0.591052172672043"
Il y a plus d’overfitting (écart de 1 point entre train et test) dans le GBM que dans le GLM, mais il y a aussi un fit de bien meilleur qualité, on va comparer les courbes GLM et GBM sur le test.
pred_gbmtest%>%
na.omit%>%
arrange(-pred)%>%
mutate(y=cumsum(obs)/sum(obs),x=(1:nrow(.))/nrow(.))%>%
{
data.frame(y100=quantile(.$y,0:100/100),
x100=quantile(.$x,0:100/100))
}->Lorenz_gbm_points
g <- g +
geom_line(data = Lorenz_gbm_points,
aes(x=x100,y=y100,color="gbm"))
g%>%ggplotly
vars=summary(gbm_model)$var%>%as.character
summary(gbm_model)
## var rel.inf
## NBMENFISC14 NBMENFISC14 34.9627563
## NBPERSMENFISC14 NBPERSMENFISC14 18.8350015
## TP60AGE614 TP60AGE614 16.1073117
## TP60AGE114 TP60AGE114 4.9345014
## TP60AGE514 TP60AGE514 2.8233126
## RD14 RD14 2.6426279
## TP60TOL114 TP60TOL114 1.6481053
## PPLOGT14 PPLOGT14 1.5172798
## PIMPOT14 PIMPOT14 1.5128370
## TP60AGE414 TP60AGE414 1.4323254
## PBEN14 PBEN14 1.2352326
## D114 D114 1.2160091
## TP60TOL214 TP60TOL214 1.1362209
## PPAT14 PPAT14 1.1027992
## TP60AGE214 TP60AGE214 1.1016187
## PPEN14 PPEN14 1.0699677
## TP60AGE314 TP60AGE314 1.0027726
## D914 D914 0.8113241
## PRA14 PRA14 0.8075889
## PTSAC14 PTSAC14 0.7694872
## PPFAM14 PPFAM14 0.6743689
## PPMINI14 PPMINI14 0.6011457
## TP6014 TP6014 0.6000001
## PIMP14 PIMP14 0.5668264
## MED14 MED14 0.4791458
## PPSOC14 PPSOC14 0.4094331
Attention, contrairement à un GLM, ces courbes ne s’interprètent pas comme des effets toutes choses égales par ailleurs parce qu’elles s’appuient sur la distribution réelle des données.
plot(gbm_model,i.var=vars[1])
plot(gbm_model,i.var=vars[2])
plot(gbm_model,i.var=vars[3])
plot(gbm_model,i.var=vars[4])
A comparer avec les valeurs réellement prises par ces variables, les extrêmes sont estimés mais souvent les valeurs sont aberrantes et non significatives.
## Monotonie des variables
Si nécessaire on peut imposer la monotonie des variables, ceci réduira le surapprentissage et donnera des effets plus propres/plus interprétables.
dans gbm : var.monotone il convient de donner un vecteur nommé (avec les noms des variables explicatives) et des valeurs -1 si décroissante, 1 si croissante, 0 si absence de contrainte.
monotony=rep(0,ncol(data_model)-1)
names(monotony) <- setdiff(names(data_model),"nb_ET")
monotony["TP60AGE114"] <- 1
params=c(shrinkage=.01,
nb_trees=1000,
depth=10,
nminobs=10,
bag.frac=.2)
gbm_model <- gbm(nb_ET~.
,var.monotone = monotony
,data=data_model[train,]
,verbose=T
,train.fraction=.8
,shrinkage = params[1]
,n.trees = params[2]
,interaction.depth = params[3]
,n.minobsinnode = params[4]
,bag.fraction = params[5]
)
## Distribution not specified, assuming gaussian ...
## Iter TrainDeviance ValidDeviance StepSize Improve
## 1 161.4961 125.0582 0.0100 1.2216
## 2 159.8875 123.5936 0.0100 1.5888
## 3 158.4206 122.4782 0.0100 1.4064
## 4 157.5669 121.8906 0.0100 0.7681
## 5 155.6164 120.4929 0.0100 1.5742
## 6 154.0022 119.2011 0.0100 1.5202
## 7 152.0423 117.4708 0.0100 1.8606
## 8 151.0418 116.6187 0.0100 1.0506
## 9 149.9695 115.5649 0.0100 1.1009
## 10 148.6748 114.5025 0.0100 1.2641
## 20 135.6161 103.5849 0.0100 1.1795
## 40 115.8605 87.1145 0.0100 0.3968
## 60 101.1442 76.1938 0.0100 0.5126
## 80 91.3902 69.6246 0.0100 0.4239
## 100 85.1180 66.8352 0.0100 0.2402
## 120 80.5801 64.0562 0.0100 -0.0309
## 140 76.0434 62.3596 0.0100 0.3002
## 160 71.9982 61.3791 0.0100 -0.0483
## 180 69.7115 61.4357 0.0100 0.0863
## 200 67.8822 61.7476 0.0100 0.0993
## 220 66.3625 61.9588 0.0100 -0.0540
## 240 64.6558 62.1391 0.0100 0.0443
## 260 63.4095 62.4444 0.0100 -0.0294
## 280 62.4027 62.7356 0.0100 -0.0323
## 300 61.4847 63.2967 0.0100 -0.0543
## 320 60.1014 64.0942 0.0100 -0.0077
## 340 59.6189 64.2372 0.0100 0.0065
## 360 59.1495 64.5275 0.0100 -0.0536
## 380 58.2259 65.3635 0.0100 -0.0273
## 400 57.5878 65.7638 0.0100 0.0018
## 420 56.9295 66.4420 0.0100 -0.0161
## 440 56.3612 66.8646 0.0100 -0.0596
## 460 55.7183 67.6218 0.0100 -0.0187
## 480 55.3136 67.7511 0.0100 -0.0235
## 500 54.9858 67.5345 0.0100 -0.0296
## 520 54.3594 67.7049 0.0100 0.0438
## 540 53.9246 67.7988 0.0100 -0.0837
## 560 53.5223 68.1030 0.0100 0.0430
## 580 53.2301 67.7800 0.0100 -0.0834
## 600 52.8462 67.7087 0.0100 -0.0289
## 620 52.3831 67.9810 0.0100 -0.0006
## 640 51.9556 68.2019 0.0100 -0.0276
## 660 51.4675 68.1175 0.0100 -0.0459
## 680 51.2878 67.7366 0.0100 0.0050
## 700 50.7660 68.0417 0.0100 -0.0207
## 720 50.4559 67.6265 0.0100 -0.0792
## 740 49.9923 68.1120 0.0100 0.0101
## 760 49.6616 68.2458 0.0100 0.0180
## 780 49.4424 68.4961 0.0100 -0.0425
## 800 49.0898 68.6180 0.0100 0.0125
## 820 48.6806 68.4983 0.0100 -0.0262
## 840 48.4303 68.7273 0.0100 -0.0407
## 860 48.0850 69.3697 0.0100 -0.0104
## 880 47.6589 69.4833 0.0100 -0.0119
## 900 47.3868 69.9136 0.0100 -0.0164
## 920 47.0484 70.5000 0.0100 -0.0194
## 940 46.6978 70.2430 0.0100 -0.0012
## 960 46.5355 70.8139 0.0100 -0.0087
## 980 46.1201 70.6895 0.0100 -0.0211
## 1000 45.6229 70.8181 0.0100 -0.0228
pred_gbmtrain=data.frame(pred=predict(gbm_model,
newdata = data_model[train,]),
obs=target_train)
## Using 166 trees...
pred_gbmtest=data.frame(pred=predict(gbm_model,
newdata = data_model[-train,]),
obs=target_test)
## Using 166 trees...
print(paste("gbm test:",Gini(pred_gbmtest)))
## [1] "gbm test: 0.845622319906236"
print(paste("gbm train:",Gini(pred_gbmtrain)))
## [1] "gbm train: 0.861493418882915"
print(paste("glm:",Gini(pred_glmnettest)))
## [1] "glm: 0.591052172672043"
plot(gbm_model,i.var="TP60AGE114")
Pour en savoir plus sur la simplification des modèles de machine learning complexes, voici un article sur “interpretable machine learning” IML
Dans le cas spécifique des GBM des méthodes adhoc existent parce qu’un arbre peut être vu comme un ensemble d’indicatrices sur des croisements de variables.
Il y a potentiellement autant de croisements que la profondeur de l’arbre.
Si on décide d’utiliser des “stumps” arbres de profondeur 1 alors chaque arbre est une indicatrice sur 1 variable. Autrement dit ce sont des GLM (parce qu’il y a quand même une fonction de lien) avec des interactions d’ordre 1.
Si on décide d’utiliser des arbres de profondeur 2 alors les arbres construisent des indicatrices sur des croisements de 2 variables. Autrement dit ce sont des GLM avec des interactions d’ordre 2.
Or un GBM (idem pr random forest) est une somme pondérée d’arbres de décisions.
Et une somme pondérée (combinaison linéaire) de GLM est encore un GLM.
Bref on peut ramener un GBM de profondeur 1 ou 2 à un GLM avec des interactions d’ordre 1 ou 2.
Ceci dit il y aura beaucoup de coefficients donc on est plus proche d’une GAM non paramétrique que d’un GLM (au lieu de splines on utilise des fonctions constantes par morceaux)